Continue working with the dataset on the effects of global warming on hibernation

Iris Pretzlaff at Hamburg University studied hibernation of dormice with the goal of understanding how global warming will affect physiology during hibernation and how this can affect population viability

Hypotheses:

The effect of temperature on metabolic rate differs among individuals in torpor and in normothermic state
- A quadratic effect of temperature on metabolic rate of torpid individuals is expected, so that the metabolic rate is lowest during optimal temperature range and increases for too low / too high temperatures
- No effect of temperature on metabolic rate of normothermic individuals is expected because they should be able to thermoregulate
Probability of arousal during hibernation is related to ambient temperature and time since hibernation. For both variables a quadratic effect is expected, so that the probability of arousal is lower at the intermediate values of temperature and time since arousal and increases for low / high temperature and the beginning / end of hibernation
Body mass is expected to decrease with high arousal frequency and high average metabolic rate one month prior to body mass measurement

Data collection

Iris observed 12 and 13 individuals over winters of 2009 and 2011 in wooden boxes, kept outside in their natural conditions. Skin temperature was measured with data loggers. Torpor and arousal bouts were assigned based on skin temperature: Torpor was defined as any time when skin temperature decreased < 30°C (Tskin of active, normothermic animals is around 36°C) and an arousal when skin temperature increased > 30°C between torpor bouts, or when the difference Tskin-Tnest was > 30°C.

Metabolic rate was measured as oxygen consumption (ml \(O_{2}\)/h) using open-flow respirometry. For more details see Pretzlaff et al. (2021).

Focus: Probability of arousal

We focus on testing the hypothesis about the probability of arousal.

Exploratory plots

xyplot(state ~ Tamb|ID, data = df)

xyplot(state ~ HoursBegin|ID, data = df)

s2009 <- df[df$Year %in% c(2009, 2010), ]
s2011 <- df[df$Year %in% c(2011, 2012), ]

xyplot(state ~ idate | ID, data = s2009)

xyplot(state ~ idate | ID, data = s2011)

More data preparation

We rescale the data for better model convergence. Feel free to try and fit the model without this step (I already did that) and see what happens.

df <- df[! (is.na(df$HoursBegin) | is.na(df$Tamb)), ]
min(df$HoursBegin)

## [1] 19

# coding year as a factor
df$Y <- 2009
df$Y[df$Year %in% c(2011, 2012)] <- 2011
df$Y <- as.factor(df$Y)

### rescaling time since hibernation and temperature for better convergence
df$HoursScale <- (df$HoursBegin - mean(df$HoursBegin, na.rm = T)) / sd(df$HoursBegin, na.rm = T)
df$HoursScale2 <- df$HoursScale^2
hist(df$HoursScale, 
     xlab = 'Standarized time since begin hibernation',
     main = NULL, col = 'black')

df$T_scale <- (df$Tamb - mean(df$Tamb, na.rm = T)) / sd(df$Tamb, na.rm = T)
df$Tscale2 <- df$T_scale ^ 2
hist(df$T_scale, xlab = 'Standardized temperature', 
     main = NULL, col ='black')

Fit an unrealistic model, a GLM with Binomial distribution and logit link.

mod0_probA <- glm(state ~ HoursScale + HoursScale2 + 
                    T_scale + Tscale2 + Y, 
                  family = binomial(link = 'logit'),
                  data = df)
summary(mod0_probA)

## 
## Call:
## glm(formula = state ~ HoursScale + HoursScale2 + T_scale + Tscale2 + 
##     Y, family = binomial(link = "logit"), data = df)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.9353  -0.5560  -0.3860  -0.2337   2.7603  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -3.35664    0.05807 -57.803  < 2e-16 ***
## HoursScale   0.01895    0.02040   0.929    0.353    
## HoursScale2  0.72028    0.02783  25.883  < 2e-16 ***
## T_scale      0.45224    0.03498  12.930  < 2e-16 ***
## Tscale2      0.11389    0.01913   5.954 2.61e-09 ***
## Y2011        0.95357    0.04988  19.119  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 14940  on 16579  degrees of freedom
## Residual deviance: 11951  on 16574  degrees of freedom
## AIC: 11963
## 
## Number of Fisher Scoring iterations: 6

Any idea how we can conduct model diagnostics for this model? What about calculating dispersion, do you remember how we can do a quick check whether overdispersion is potentially an issue?

Fit a mixed-effect model with random intercept.

Here we assume that time since the beginning of hibernation and temperature have same effects across all individuals, but that there is some random variation among individuals in the intercept for i) time since begin of hibernation and ii) temperature for each individual. We use glmmTMB package to fit these models because we want to be able to plot the predictions in the end (+ ideally fit autocorrelation per individual)

modRandID_prob <- glmmTMB(state ~ HoursScale + HoursScale2 +
                          T_scale + Tscale2 +
                          Y + (1|ID), 
                         family = binomial,
                         data=df)
summary(modRandID_prob)

##  Family: binomial  ( logit )
## Formula:          
## state ~ HoursScale + HoursScale2 + T_scale + Tscale2 + Y + (1 |      ID)
## Data: df
## 
##      AIC      BIC   logLik deviance df.resid 
##  11147.4  11201.4  -5566.7  11133.4    16573 
## 
## Random effects:
## 
## Conditional model:
##  Groups Name        Variance Std.Dev.
##  ID     (Intercept) 0.6676   0.8171  
## Number of obs: 16580, groups:  ID, 25
## 
## Conditional model:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -3.55700    0.23680 -15.021  < 2e-16 ***
## HoursScale   0.11709    0.02278   5.141 2.74e-07 ***
## HoursScale2  0.80785    0.03111  25.969  < 2e-16 ***
## T_scale      0.48802    0.03645  13.387  < 2e-16 ***
## Tscale2      0.14762    0.01906   7.746 9.49e-15 ***
## Y2011        0.92287    0.33263   2.774  0.00553 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Diagnostics of the fitted random-intercept model. Let us use DHARMa package.

sim <- simulateResiduals(modRandID_prob, plot = T)

testUniformity(sim, plot = F)

## 
##  Asymptotic one-sample Kolmogorov-Smirnov test
## 
## data:  simulationOutput$scaledResiduals
## D = 0.011804, p-value = 0.0197
## alternative hypothesis: two-sided

testQuantiles(sim, plot = T)

## 
##  Test for location of quantiles via qgam
## 
## data:  simulationOutput
## p-value < 2.2e-16
## alternative hypothesis: both

testDispersion(sim)

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 0.92247, p-value = 0.296
## alternative hypothesis: two.sided

# plot residuals vs our predictors
plotResiduals(sim, form = df$T_scale)

plotResiduals(sim, form = df$HoursScale)

What does this diagnostics plots tell us?

Test for temporal autocorrelation since the measurements were taken over time.

sim_ID10 <- recalculateResiduals(sim, sel = df$ID == 10)
testTemporalAutocorrelation(sim_ID10, time = df$Timecon[df$ID == '10'])

## 
##  Durbin-Watson test
## 
## data:  simulationOutput$scaledResiduals ~ 1
## DW = 1.6427, p-value = 1.239e-07
## alternative hypothesis: true autocorrelation is not 0

Let us use a similar function we used in the lecture to test temporal autocorrelation for each individual in the dataset.

levels(df$ID)

##  [1] "10"   "11A"  "11B"  "12"   "14A"  "14B"  "2"    "21"   "22"   "2222"
## [11] "23"   "24A"  "24B"  "38A"  "38B"  "38C"  "4"    "57"   "64"   "70"  
## [21] "79"   "86"   "97A"  "97B"  "98"

fun_autocor(sim, df, ID = '10')

## [[1]]
## 
##  Durbin-Watson test
## 
## data:  simulationOutput$scaledResiduals ~ 1
## DW = 1.6427, p-value = 1.239e-07
## alternative hypothesis: true autocorrelation is not 0
## 
## 
## [[2]]
##      DurbW           p ID
## DW 1.64269 1.23855e-07 10

Test autocorrelation in all studies

dat_autoc <- data.frame(DurbW =numeric(), p = numeric())
dat_autoc <- rbind(dat_autoc, fun_autocor(sim, df, 10)[[2]])
for(i in unique(df$ID)[-1]){
  dat_autoc <- rbind(dat_autoc, fun_autocor(sim, df, ID = i)[[2]])
}

dat_autoc$sign <- ifelse(dat_autoc$p < 0.05,'Yes','No')
sum(dat_autoc$sign == 'Yes')

## [1] 25

length(levels(df$ID))

## [1] 25

So, autocorrelation is an issue

Fit the Binomial model with autocorrelation

Attention: we assume we have time running continuously for each individual. But in fact, for some of them there are little gaps in the records (for temperature and metabolic rate)

df$Time_fac <- as.factor(df$Timecon)
modRandAR_prob <- glmmTMB(state ~ HoursScale + HoursScale2 +
                          T_scale + Tscale2 +
                          Y + (1|ID) +  ar1(Time_fac - 1|ID), 
                         family = binomial,
                         data=df)
summary(modRandAR_prob)

##  Family: binomial  ( logit )
## Formula:          
## state ~ HoursScale + HoursScale2 + T_scale + Tscale2 + Y + (1 |  
##     ID) + ar1(Time_fac - 1 | ID)
## Data: df
## 
##      AIC      BIC   logLik deviance df.resid 
##   3893.1   3962.6  -1937.6   3875.1    16571 
## 
## Random effects:
## 
## Conditional model:
##  Groups Name        Variance  Std.Dev.  Corr      
##  ID     (Intercept) 8.467e-07 9.201e-04           
##  ID.1   Time_fac1   1.937e+03 4.401e+01 0.93 (ar1)
## Number of obs: 16580, groups:  ID, 25
## 
## Conditional model:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -60.3817     6.3070  -9.574  < 2e-16 ***
## HoursScale    1.1673     0.9665   1.208  0.22710    
## HoursScale2  13.7559     1.6632   8.271  < 2e-16 ***
## T_scale       8.3990     1.3541   6.203 5.55e-10 ***
## Tscale2       1.6549     0.5481   3.020  0.00253 ** 
## Y2011        19.2246     3.1915   6.024 1.71e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

How do you interpret these results? What do you think about this model compared to the one without autocorrelation?

Model diagnostics for the model with autocorrelation

simAR1 <- simulateResiduals(modRandAR_prob, plot = T)

## DHARMa:testOutliers with type = binomial may have inflated Type I error rates for integer-valued distributions. To get a more exact result, it is recommended to re-run testOutliers with type = 'bootstrap'. See ?testOutliers for details

testUniformity(simAR1, plot = F)

## 
##  Asymptotic one-sample Kolmogorov-Smirnov test
## 
## data:  simulationOutput$scaledResiduals
## D = 0.066222, p-value < 2.2e-16
## alternative hypothesis: two-sided

testQuantiles(simAR1, plot = T)

## 
##  Test for location of quantiles via qgam
## 
## data:  simulationOutput
## p-value < 2.2e-16
## alternative hypothesis: both

testDispersion(simAR1)

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 0.6603, p-value < 2.2e-16
## alternative hypothesis: two.sided

# plot residuals vs our predictors
plotResiduals(simAR1, form = df$T_scale)

plotResiduals(simAR1, form = df$HoursScale)

What are main issues with model diagnostics? What can be done to fit a model that would match the data structure better?
We have fitted the models assuming the slopes are exactly the same across individuals and only intercepts differ. How realistic is that? Should we have tried a different random structure?

From the note on DHARMa documentation:

Standard DHARMa simulations from models with (temporal / spatial / phylogenetic) conditional autoregressive terms will still have the respective temporal / spatial / phylogenetic correlation in the DHARMa residuals, unless the package you are using is modelling the autoregressive terms as explicit REs and is able to simulate conditional on the fitted REs. This has two consequences: 1. If you check the residuals for such a model, they will still show significant autocorrelation, even if the model fully accounts for this structure. 2. Because the DHARMa residuals for such a model are not statistically independent any more, other tests (e.g. dispersion, uniformity) may have inflated type I error, i.e. you will have a higher likelihood of spurious residual problems.

Attention! This is very slow and takes long time to run. I suggest you run it after the session

##DO not run
simAR1_rot <- simulateResiduals(modRandAR_prob, plot = T, rotation = 'estimated')
testUniformity(simAR1_rot, plot = F)
testQuantiles(simAR1_rot, plot = T)

testDispersion(simAR1_rot)

# plot residuals vs our predictors
plotResiduals(simAR1_rot, form = df$T_scale)
plotResiduals(simAR1_rot, form = df$HoursScale)

Conclusion: know what function is doing exactly. Do not just copy-paste the codes!

Model inference

To test the hypothesis that both temperature and time have a quadratic effect on probability of arousal, we fit the models without each quadratic term in turn, and run LRTs.

modRandAR_no2T_prob <- glmmTMB(state ~ HoursScale + HoursScale2 +
                          T_scale + 
                          Y + (1 | ID) +  ar1(Time_fac - 1|ID), 
                         family = binomial,
                         data=df)

## Warning in fitTMB(TMBStruc): Model convergence problem; singular convergence
## (7). See vignette('troubleshooting')

modRandAR_no2Hour_prob <- glmmTMB(state ~ HoursScale + 
                          T_scale + Tscale2 +
                          Y + (1 | ID) +  ar1(Time_fac - 1|ID), 
                         family = binomial,
                         data=df)

## Warning in fitTMB(TMBStruc): Model convergence problem; false convergence (8).
## See vignette('troubleshooting')

anova(modRandAR_prob, modRandAR_no2T_prob)  ## significant quadratic effect

## Data: df
## Models:
## modRandAR_no2T_prob: state ~ HoursScale + HoursScale2 + T_scale + Y + (1 | ID) + ar1(Time_fac - , zi=~0, disp=~1
## modRandAR_no2T_prob:     1 | ID), zi=~0, disp=~1
## modRandAR_prob: state ~ HoursScale + HoursScale2 + T_scale + Tscale2 + Y + (1 | , zi=~0, disp=~1
## modRandAR_prob:     ID) + ar1(Time_fac - 1 | ID), zi=~0, disp=~1
##                     Df    AIC    BIC  logLik deviance Chisq Chi Df Pr(>Chisq)
## modRandAR_no2T_prob  8 3905.8 3967.5 -1944.9   3889.8                        
## modRandAR_prob       9 3893.1 3962.6 -1937.6   3875.1 14.65      1  0.0001294
##                        
## modRandAR_no2T_prob    
## modRandAR_prob      ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Note we have ‘Warning about false convergence’ for the model without the quadratic effect of time. Inspect its output and think what are the possible causes and what can be done?

summary(modRandAR_no2Hour_prob)

##  Family: binomial  ( logit )
## Formula:          
## state ~ HoursScale + T_scale + Tscale2 + Y + (1 | ID) + ar1(Time_fac -  
##     1 | ID)
## Data: df
## 
##      AIC      BIC   logLik deviance df.resid 
##   3997.9   4059.7  -1991.0   3981.9    16572 
## 
## Random effects:
## 
## Conditional model:
##  Groups Name        Variance  Std.Dev.  Corr      
##  ID     (Intercept) 1.661e-06  0.001289           
##  ID.1   Time_fac1   1.242e+03 35.242402 0.93 (ar1)
## Number of obs: 16580, groups:  ID, 25
## 
## Conditional model:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -35.4141     3.8188  -9.274  < 2e-16 ***
## HoursScale   -2.9674     0.9732  -3.049  0.00229 ** 
## T_scale      11.2666     1.2090   9.319  < 2e-16 ***
## Tscale2       2.2701     0.3733   6.081 1.19e-09 ***
## Y2011        12.4581     2.7251   4.572 4.84e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Plot predictions

Prepare a new data frame. First predict on the scale of the linear predictor.

min(df$T_scale)

## [1] -4.340506

min(df$HoursScale)

## [1] -2.060011

new_dat <- expand.grid('ID' = unique(df$ID), 
                       'Y' = unique(df$Y),
                       T_scale = seq(from = min(df$T_scale), to = 
                                    max(df$T_scale), length.out = 20),
                       HoursScale = seq(from = min(df$HoursScale), to = 
                                    max(df$HoursScale), length.out = 20))

new_dat$HoursScale2 <- new_dat$HoursScale^2
new_dat$Tscale2 <- new_dat$T_scale^2
head(new_dat)

##    ID    Y   T_scale HoursScale HoursScale2  Tscale2
## 1  10 2011 -4.340506  -2.060011    4.243644 18.83999
## 2 11A 2011 -4.340506  -2.060011    4.243644 18.83999
## 3 11B 2011 -4.340506  -2.060011    4.243644 18.83999
## 4  12 2011 -4.340506  -2.060011    4.243644 18.83999
## 5 14A 2011 -4.340506  -2.060011    4.243644 18.83999
## 6 14B 2011 -4.340506  -2.060011    4.243644 18.83999

new_dat$Time_fac <- as.factor(rep('1', nrow(new_dat)))  ## keeping just to one time point, to avoid variaiton here
new_dat$pred_fixTMB <- predict(modRandAR_prob, newdata = new_dat, re.form = NA, type = 'link')

new_dat <- new_dat[order(new_dat$ID, new_dat$Y, new_dat$T_scale), ]
ggplot(data = new_dat, aes(x = T_scale, y = pred_fixTMB, group = as.factor(HoursScale), col = as.factor(HoursScale))) +
  geom_point() + theme_bw() +
  facet_grid(cols = vars(Y)) +
  xlab('Standardized temperature') + ylab('Linear predictor')

ggplot(data = new_dat, aes(x = HoursScale, y = pred_fixTMB, group = as.factor(HoursScale), col = as.factor(T_scale))) +
  geom_point() + theme_bw() +
  facet_grid(cols = vars(Y)) +
  xlab('Standardized time sing begin of hibernation') + ylab('Linear predictor')

What is the probability to arouse close to the end of hibernation

max(df$HoursScale)

## [1] 1.841177

new_dat <- expand.grid('ID' = unique(df$ID), 
                       'Y' = unique(df$Y),
                       T_scale = seq(from = min(df$T_scale), to = 
                                    max(df$T_scale), length.out = 20),
                       HoursScale = 1.6)  ## set time to median, to avoid additional variaiton

new_dat$HoursScale2 <- new_dat$HoursScale^2
new_dat$Tscale2 <- new_dat$T_scale^2
head(new_dat)

##    ID    Y   T_scale HoursScale HoursScale2  Tscale2
## 1  10 2011 -4.340506        1.6        2.56 18.83999
## 2 11A 2011 -4.340506        1.6        2.56 18.83999
## 3 11B 2011 -4.340506        1.6        2.56 18.83999
## 4  12 2011 -4.340506        1.6        2.56 18.83999
## 5 14A 2011 -4.340506        1.6        2.56 18.83999
## 6 14B 2011 -4.340506        1.6        2.56 18.83999

new_dat$Time_fac <- as.factor(rep('1', nrow(new_dat)))  ## keeping just to one time point, to avoid variaiton here
new_dat$pred_fixTMB <- predict(modRandAR_prob, newdata = new_dat, re.form = NA, type = 'response')

ggplot(data = new_dat, aes(x = T_scale, y = pred_fixTMB)) +
  geom_line() + theme_bw() +
  facet_grid(cols = vars(Y)) +
  xlab('Standardized temperature') + ylab('Probability of arousal')

Add CI to the prediction on probability to arouse from hibernation

new_dat$pred_fixTMB <- predict(modRandAR_prob, newdata = new_dat, re.form = NA, type = 'response')
new_dat$pre_SE_fixTMB <- predict(modRandAR_prob, newdata = new_dat, re.form = NA, se.fit = TRUE, type = 'response')$se.fit 
new_dat$pred_upper_fix <- new_dat$pred_fixTMB + qnorm(0.975)*new_dat$pre_SE_fixTMB 
new_dat$pred_lower_fix <- new_dat$pred_fixTMB - qnorm(0.975)*new_dat$pre_SE_fixTMB 
head(new_dat)

##    ID    Y   T_scale HoursScale HoursScale2  Tscale2 Time_fac  pred_fixTMB
## 1  10 2011 -4.340506        1.6        2.56 18.83999        1 8.682929e-05
## 2 11A 2011 -4.340506        1.6        2.56 18.83999        1 8.682929e-05
## 3 11B 2011 -4.340506        1.6        2.56 18.83999        1 8.682929e-05
## 4  12 2011 -4.340506        1.6        2.56 18.83999        1 8.682929e-05
## 5 14A 2011 -4.340506        1.6        2.56 18.83999        1 8.682929e-05
## 6 14B 2011 -4.340506        1.6        2.56 18.83999        1 8.682929e-05
##   pre_SE_fixTMB pred_upper_fix pred_lower_fix
## 1  0.0007073774    0.001473263   -0.001299605
## 2  0.0007073774    0.001473263   -0.001299605
## 3  0.0007073774    0.001473263   -0.001299605
## 4  0.0007073774    0.001473263   -0.001299605
## 5  0.0007073774    0.001473263   -0.001299605
## 6  0.0007073774    0.001473263   -0.001299605

# predictions for the mean response across species
ggplot(data = new_dat, aes(x = T_scale, y = pred_fixTMB)) +
  theme_bw() + facet_grid(cols = vars(Y)) +
  geom_ribbon(data = new_dat, aes(x = T_scale, ymin = pred_lower_fix, ymax = pred_upper_fix), alpha = 0.1) +
  geom_line(lty  = 3) +
    xlab('Standardized temperature') + ylab('Probability of arousal')

Attention: we only used predictions with the fixed part of the model, not including random effects at all!

Repetition: how do you go from the scale of the response to the linear predictor?

new_dat$pred_fixTMB_link <- predict(modRandAR_prob, newdata = new_dat, re.form = NA, type = 'link')
new_dat$pred_trans <- exp(new_dat$pred_fixTMB_link)/(1 + exp(new_dat$pred_fixTMB_link))

ggplot(data = new_dat, aes(x = T_scale, y = pred_trans)) +
  geom_line() + theme_bw() +
  facet_grid(cols = vars(Y)) +
  xlab('Standardized temperature') + ylab('Probability of arousal')

Practical session: Generalised Linear Mixed-Effects models

Radchuk Viktoriia

2022-1-22